C.................... RSMNSD.FOR ..................................... 
C      PROGRAM FLIP
C.... this main module for the Field Line Interhemispheric Plasma (FLIP)
C.... model. It reads initial conditions and sets up the
C.... coordinate system and neutral atmosphere, note s(4500) for d300
C
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER ION,NEQ,WUNIT8,IRIHED,ICAUTN,IVERT,IVPERP,NCOL_EXB,JDAY
     > ,UTORLT,NOCON,AUR_FILE,APFILE,JZLO,JJ,FSTART,JZ,JBLATPR,J250,I
     > ,FEXISTS 
      INTEGER NUMUNIT,UNITNUMS(22)    !.. Unit numbers for writing rate warning
      REAL SEC,DEC,ETRAN,F107,F107A,AP,BLON,D,T,CHI,SATN,SATF,SAT,EUV
     > ,UVFAC,GLON,GLONN,GLONF,CHIN,CHIF,GLATN,GLATF,GLATJ,AP3HR(33001)
     > ,GRD,RATFAC,SW(25),UTHRS,UHED(2),SATEQ,L_STAG,L_SHORT
     > ,RE,PLAT,PLON,PI,ANGVEL,TIMEON,TIMEOF,AURALT,ZTOPTE,L_PHASE,SX
     > ,TD_PHASE,DELTIM,PHASE,SECWF8,GLOND,GLATD,BCOMPU,KPREAL,ZFLUX(37)
     > ,SATX,CHIX,GLONX,GLATX,SATNSAV
      !... TN,OX,O2,N2,HE,H perturbations from GWs, Nth, East, Ver, Bmer winds
      REAL GWDENTN(9),GWIND(3),GW_BMER   !.. GW added 2023-02-06

     
      REAL UMZ(2,2),AP_TO_KP(401)
      REAL EFAC,EFACN,EFACS,ECDT,ECHROM,ECORON   !.. Solar eclipse factors and time step
      REAL TWISAV                                !.. Save twilight time step
      DOUBLE PRECISION FYLIM,HPLOSS,ZTRAN        !.. For calculating limiting flux
      DIMENSION S(50000),D(19),T(2),FD(9),COSDIP(401),Q(401)
     >   ,V(2),COOL(22)
      COMMON/V91/RVBT(401)
      COMMON/PALT/ZLO,ZHI,JWR,JTI
      COMMON/RFAC/RATFAC(99)
      COMMON/RK/VC(401)
      !----- common to hold the total ion production rates
      REAL SUMION,SUMEXC
      COMMON/BPRTOT/SUMION(3,12,401),SUMEXC(3,12,401)
      !.. N= [O+], [H+], sum minor, [He+] densities, Ti= H+,O+,e temps 
      !.. F= cty eqn functions
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      !.. Ap, solar declination, ephemeris transit, magnetic longitude,
      !.. solar indices, UT in secs, magnetic latitude
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      !.. O+,H+ velocities, field factor (obsolete), magnetic field, 
      !.. GR=gravity, R=geocentric distance,SL= distance along B
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      !.. [O],[H],[N2],[O2], O+ photoionization rate, Tn
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      !.. used in minor ions [O+], V(O+), [He], ???, He+ prod, ???
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      !.. Altitude(km), time step, ???, ???, ???, ITF= change solution
      !.. mode, JMAX=# grid points, JMAX1=JMAX-1,# Newton iterations, # ions
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      !.. NSAVE= saved [O+] and [H+], O+, H+ flux, EHT=thermal e heating rate
      COMMON/SAV/NSAVE(2,401),TISAV(3,401),FY(2,401),UN(401),EHT(3,401)
      !.. solar zenith angle, ??? ??? ??? ???, max # time steps ??? ???
      COMMON/LPS/SZA(401),EPSN,DC,ISPC,I1,I2,IPV,JTIMAX,IPP,ISKP
      !.. Stores vibrational N2 and odd N densities for printing. 
      !.. RCON= effect of vibrational N2 on O+ + N2 reaction
      COMMON/STAW/STR(12,401),RCON(401)
      !.. Factors for multiplying EUV, UV, and MSIS
      COMMON/SOL/UVFAC(59),EUV
      !.. RE= Earth radius, (PLAT,PLON) = geographic location of north
      !.. magnetic pole, Pi, angular velocity of Earth 
      COMMON/GPAR/RE,PLAT,PLON,PI,ANGVEL
      !.. Minor ion and neutral densities O+(2D) O+(2P) N(2D) N(4S), NO 
      COMMON/MIB/OP2D(401),OP2P(401),N2D(401),N4S(401),NNO(401)
      !.. densities, velocities, and masses for minor ions solved by Newton
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)
      !.. Common for O+ - O collision frequencies
      COMMON/CFACT/DHDU_FAC,COLFAC,CFAC(9),ISCH(9)
      !.. Switches for wind type and equations solved
      COMMON/CONFAC/IWIND,IVIBN2,IODDN,INPLUS,IHEP,ITEMP,IHPOP
      !.. minor neutral species on field line grid
      COMMON/VODDN/VN2D(401),VN4S(401),VNNO(401),VO1D(401),IVERT

      !.. JB= grid point of lower boundary, 3 hour Aps, winds flag, SW=MSIS
      !.. switches, ZPROD = upper boundary for EUV ion production
      DATA JB/1/, AP3HR/-1,33000*0/, IWTIM/-1/,SW/15*1,10*1/,ZPROD/700./
      !.. ITSAV ITAU track elapsed time, ISOL= type of solution, HFAC=
      !.. e heating factor, MAXIT= max Newton iters, JTIDT used to set time
      !.. ZLBDY=lower boundary altitude for solving equations, IRIHED IRI/HWM
      !.. wind switch, max solar zenith angle for ion production
      DATA ITSAV,ITAU, ISOL, HFAC, MAXIT, JTIDT, ZLBDY,  IRIHED, SETCHI
     >  / -10000, 0  ,  0  ,  0.0 ,  17  ,   5   , 120. ,  0   ,  1.7/
      !.. north/south lat and long 250 km above station, Z0=absolute lower
      !.. boudary altitude, SCAL, GRD scale grid B and vertical grid points
      DATA GLONN, GLONF, GLATN, GLATF, Z0, SCAL, GRD / 7*0.0/,APFILE/-1/
      DATA JHFLX/0/,WUNIT8/1/,NCOL_EXB/0/, JJJ/0/, ICAUTN/0/
      DATA IVPERP/-1/,SCALK/-1.0/,PCO_MAX/0.0/, ISYM/0/, SECWF8/-999/
      DATA HMF2N,HMF2S/250.0,250.0/, JK,JC/10,10/, NOCON/0/, JBLATPR/0/
      DATA PCO_LIM/1.03/   !.. minimum allowable L-value
      DATA PCO_WIND/1.08/  !.. minimum allowable L-value for vertical grid
      DATA EFAC,EFACN,EFACS,ECDT/4*1.0/   !.. Eclipse factors
      !.. Unit numbers for writing reaction rate warning
      DATA NUMUNIT/9/,UNITNUMS/3,4,6,9,12,13,16,56,66,13*0/
      DATA FEXISTS/-9/    !.. Used in routine TFILE to detect a file
      !.. Conversion factors from Ap to Kp and vice versa.
      !..Kp = 0   0+   1-   1   1+   2-   2   2+   3-   3   3+   4-   4   4+
      !..ap = 0   2    3    4   5    6    7   9    12   15  18   22   27  32
      !..Kp = 5-   5  5+  6-   6  6+   7-   7   7+    8-   8    8+   9-   9      
      !..ap = 39  48  56  67  80  94  111  132  154  179  207  236  300  400
      DATA AP_TO_KP/1*0,1*0,1*0.3,1*0.7,1*1,1*1.3,1*1.7,2*2,3*2.3,3*2.7,
     >  3*3,4*3.3,5*3.7,5*4,7*4.3,9*4.7,8*5,11*5.3,13*5.7,14*6,17*6.3,
     >  21*6.7,22*7,25*7.3,28*7.7,29*8,64*8.3,100*8.7,1*9/
      DATA SATNSAV/24.0/
      PCO=0.0
      JC=1
      RE=6370  ! Apex Re 6378
      ANGVEL=5.2885E-9
      PI=3.14159
      IWW=0
      CALL FLIPIN(S)      !......... initialize arrays

      !.. $$$$$ UPDATE Version # 2023-01-09 here and in PHEAD & PHEADR $$$$$$
      WRITE(6,676)
 676  FORMAT(/6X,'..... FLIP MODEL with APEX Coords: ...2023-10-13 ...'/
     > /,6X,'*** This version takes 2-3 minutes to set up the Apex',
     > /,6X,' coordinates on a basic PC ***'
     > /,6X,'*** Check bottom of *LOG*.txt file for FLIP updates ***'
     > //,6X,'Consult either Phil Richards (pgrichds@gmail.com)'
     > /,6X,'or Pat Dandenault (Patrick.Dandenault@jhuapl.edu) for help'
     > /,6X,'and the appropriate level of recognition for model use.'
     > //,6X,'A FLIP run has 2 phases, The initial FLIP run gives small'
     > /,6X,'summary files of key parameters such as hmF2 and NmF2 and',
     > /,6X,'one large file. By default these output files have names'
     > /,6X,'like FR*.txt. In the 2nd PRINT phase the large FR*.txt is'
     > /,6X,'read by FLIP to produce several comprehensive output',
     > /,6X,'files with altitude profiles. By default, these files',
     > /,6X,'have names like PR*.txt.'//)
      
      CALL OPEN_FILE() !.. added by Xiang Zhao for PC

      !.. KMON= switch to subroutine for printing. ERR flag keeps
      !.. reading the headers until a proper value is obtained. If
      !.. KMON=-1, no need for large binary file (set NOWR8=1)
      NOWR8=0
 137  READ(5,*,ERR=137) KMON
      IF(KMON.EQ.-1) THEN
        NOWR8=1
        KMON=0
      ENDIF
      IF(KMON.NE.0) CALL MONPRN   !.. branch to printing routine

      WRITE(6,173)
 173  FORMAT(/6X,'--------------- BEGIN FLIP RUN ------------------'/)
      
      !.... Check if there is a file with initial densities and temperatures .......
      CALL TFILE(99,FSTART)  !.. See if a file exists
      IF(FSTART.NE.0) THEN
        ISOL=-1   !.. Signals continuation run
        CALL READ_BQ(SCALK,Q)
        CALL RRSTRT(99,PCO,Z0,SCAL,GRD,GLONN,GLONF,GLATN,GLATF,ISOL,UNN
     >  ,UNC,HMF2N,HMF2S,NMF2N,NMF2S,PLAT,PLON)
      ENDIF

      !....... Input parameters to control run. BLATPR 2020-08-30
      CALL DATRD1(ISOL,ITF,ITFS,I1,I2,IPP,MAXIT,JTIMAX,DTSET,DTMAX
     > ,JTIDT,TWIDT,ZLBDY,IREAD,HPEQ,FPAS,HFAC,IVLPS,TSTOP,JWR,ZLO,ZHI,
     > IWIND,IVIBN2,IODDN,INPLUS,IHEP,HEPRAT,ITEMP,IHPOP,IHPKP,BLATPR)
      IF(ITF.EQ.-9) IWIND=0     !.. use Hedin's wind in quick run
      TWISAV=TWIDT

      CALL RATCHK(NUMUNIT,UNITNUMS) !... check reaction rates standard
             
       !.... input new geophysical parameters
       CALL DATRD2(ISOL,IWIND,AP,BLON,F107,F107A,IDAY,SEC,UVFAC)

      !...... input parameters for setting up a new field line
      CALL DATRD3(ISOL,BLON,SEC,DTSET,F107,JMAX,PCO,Z0,SCAL,GRD
     >     ,L_STAG,L_SHORT,L_PHASE,IDAY)
      JK=20       !.. initial hmF2 point in north
      JC=JMAX-20  !.. initial hmF2 point in north

      SECSTART=SEC

       !... make sure run time is reasonable before reading F10.7
       IF(TSTOP+SEC/3600.0.GT.22*8760) THEN
         TSTOP=8785-SEC/3600
         WRITE(6,792) TSTOP
 792     FORMAT(3X,' ??* WARNING: Total run time must be <= 22 years,'
     >     ,1X,'TSTOP reset to',F10.1,' hours')
       ENDIF

      !.. Not possible to get wind from hmF2 at low or high latitudes
      !.. Need to see if files exist and if they have winds, not hmF2
      !.. 2020-07-23
      IF(PCO.LT.1.10.OR.PCO.GT.13) THEN
         IF(IABS(IWIND).GE.1) CALL TEST_WIND_FILES(IWIND,21,22)
      ENDIF

      !..... reading recently added parameters
      CALL DATRD4(ISOL,TIMEON,TIMEOF,AURALT)

      CLOSE(15)!-- close the file of run control parameters - FLIPRUN.DDD

      !...  For user input of Ap and F10.7. Convert Ap to KP for printing.
      KPREAL=0
      IF(AP(1).GT.400) AP(1)=400
      IF(AP(1).GE.0) KPREAL=AP_TO_KP(NINT(AP(1)+1))

      !... If AP(1) is negative calculate other AP from file for MSIS
      AP(2)=AP(1)
      IF(AP(1).GT.0) AP(7)=-1.0  !.. Flag indicating constant AP, F10.7
      IF(AP(1).LT.0) THEN
        CALL GET_F107AP(IDAY,SEC/3600.0,REAL(TSTOP),AP,SW,F107,F107A,
     >     KPREAL)
        APFILE=1 !.. turn on Ap from file
      ENDIF

      !.. Recalculate UVFAC and flux for variable Ap
      IF(NINT(UVFAC(58)).EQ.-1.OR.NINT(UVFAC(58)).EQ.-3)
     >  CALL FACEUV(F107,F107A,UVFAC,ZFLUX)
      IF(NINT(UVFAC(58)).EQ.-2.OR.NINT(UVFAC(58)).EQ.-4)
     >  CALL FACSEE(F107,F107A,UVFAC,ZFLUX)

      IF(ISOL.EQ.-1) ICAUTN=1
      IF(ISOL.EQ.-1) WRITE(6,820)
 820  FORMAT(/9X,' ...... Continuation of a previous FLIP run ......'/)
      IF(ISOL.EQ.0) WRITE(6,821)
 821  FORMAT(/2X,' *** Creating a new field line. Allow a few hours to'
     > ,' settle down ***')

      SECSAV=SEC
      TIMSAV=SEC
      DT=DTSET
      JMAX1=JMAX-1
      NEQ=2*JMAX-4
      IF(Z0.GE.180.AND.IVLPS.GT.0) IVLPS=0
      JQ=(JMAX+1)/2

       !..... Get solar declination and ephemeris transit time
       CALL SOLDEC(IDAY,SEC/3600.0,DEC,ETRAN)

      !..Initialize N2* cooling rate factors
      DO J=1,JMAX
        RVBT(J)=1.0
      ENDDO

      !... See if there is any ExB drift. If so, no oddN nor N2*
      CALL AMBS(1,1.3D4,0.0D0,D,T,CHI,SATEQ,GLON,GLATJ)
      IF(NINT(L_STAG).NE.0) THEN
        !.. Get maximum L at dusk for setting up proper altitude grid
        CALL GET_VPERP(IVPERP,DT,SEC/3600.,SATEQ,TSTOP,PCO,BLON,L_STAG
     >     ,L_SHORT,L_PHASE,VPERPN,NCOL_EXB,UTORLT,PCO_MIN,PCO_MAX,
     >     PCO_LIM)
         !.. warning for low L values
         IF(PCO_MIN.LT.PCO_LIM)  THEN
           WRITE(6,371) PCO_MIN,PCO_LIM
 371       FORMAT(3X,'*** WARNING: ExB drift minimum L-value =',F6.3
     >      ,1X,'is too small. Minimum L-value is reset to',F6.3/)
         ENDIF
         !.. warning for low L values and winds from hmF2
        IF(PCO_MIN.LT.PCO_WIND) THEN
          IF(IWIND.NE.0) WRITE(6,968) PCO_WIND
 968      FORMAT(5X,'CANNOT get wind from hmF2. ExB drift minimum' 
     >    ,1X,'L-value is < ',F6.3,'  using HWM-14 instead'/)
          IWIND=0.0
        ENDIF
      ENDIF

      IF(ISOL.EQ.0) CALL SET_HPEQ(PCO,HPEQ)     !.. reset the equatorial [H+]

      !.. Set up the grid points in the orthogonal q coordinates. If ExB 
      !.. drift, the field line is set up for a higher PCO to ensure
      !.. lower boundary does not rise too much
      PCO_start=PCO
      IF(PCO_MAX.GT.1) PCO_start=PCO_MAX
      IF(SCALK.LT.0) CALL QFIELD(JMAX,PCO_start,RE,Z0,SCAL,SCALK,Q)
      CALL WRITE_Q(1,JMAX,SCALK,Q)              !.. put Qs in P*.dat file
      

      !...............  MAIN TIME STEP BEGINS ....................**

      DO JTI=1,JTIMAX+1
        TWIDT=TWISAV
        !.. set the new time step  ** ++ check dL/dt in VPERP
        CALL SETDT(DT,DTSET,DTMAX,JTIDT,ITAU,Z,JMAX,ITF,ITER,TWIDT)
        IF(NINT(DTSET).EQ.59) DT=DTMAX   !..forces constant time steps from start

        !.. Calculate solar eclipse factor in North then South and reset DT
        CALL ECFUN(SEC,IDAY,45.0,ECDT,ECHROM,ECORON)
        EFACN=ECHROM  !.. 2017-09-08: Use chromospheric emission for photoionization
        IF(EFACN.LT.1) WRITE(6,'(7X,2A,F6.3,A,F6.3)') 
     >  'Eclipse in Northern Hemisphere.',
     >  '  Chromosphere irradiance fraction =',ECHROM,
     >  ';  Corona irradiance fraction =',ECORON
        CALL ECFUN(SEC,IDAY,-45.0,ECDT,ECHROM,ECORON)
        EFACS=ECHROM  !.. 2017-09-08: Use chromospheric emission for photoionization
        IF(EFACS.LT.1) WRITE(6,'(7X,2A,F6.3,A,F6.3)') 
     >  'Eclipse in Southern Hemisphere.',
     >  '  Chromosphere irradiance fraction =',ECHROM,
     >  ';  Corona irradiance fraction =',ECORON
        IF(ECDT.LT.DT.AND.(EFACN.LT.1.0.OR.EFACS.LT.1.0)) THEN
          DT=ECDT
          IF(TWIDT.GT.DT) TWIDT=DT  !.. for printinf data each time step for eclipse
        ENDIF

        !--  time to stop? Check # iterations and elapsed time
        IF(JTI.GE.JTIMAX+1) THEN
          WRITE(6,*) ' ** Completed # of time steps'
          CALL FLIP_MODS
          WRITE(6,676)
          STOP 1
        ELSEIF((SEC-SECSTART)/3600.0.GE.TSTOP) THEN
          WRITE(6,*) ' ** Completed the allotted run time'
          CALL FLIP_MODS
          WRITE(6,676)
          STOP !2
        ENDIF

      CALL ACTUAL_DAY(IDAY,SEC,JDAY,SX)    !.. gets actual day of year

        !.... Get the next AP from the file
        IF(APFILE.GE.0)
     >    CALL GET_F107AP(IDAY,SEC/3600.0,REAL(TSTOP),AP,SW,F107,F107A,
     >     KPREAL)
ccc      WRITE(6,'(22F10.1)') F107,F107A,AP

        !... adjust equatorial density for large Kp, only if from file
        IF(APFILE.GT.0.AND.IHPKP.EQ.1)
     >     CALL HP_KP(KPREAL,PCO,SEC/3600.0,JMAX,Z,N)


        !... Get ExB drift velocities. AMBS called for LT at equator
        CALL AMBS(1,1.3D4,0.0D0,D,T,CHI,SATEQ,GLON,GLATJ)
        IF(NINT(L_STAG).NE.0) THEN
          !.. determine the phase of teardrop based on Kp/Ap
          !.. L_PHASE=TD_PHASE(SEC/3600.0,AP(2),AP(3))
         CALL GET_VPERP(IVPERP,DT,SEC/3600.,SATEQ,TSTOP,PCO,BLON
     >     ,L_STAG,L_SHORT,L_PHASE,VPERPN,NCOL_EXB,UTORLT,PCO_MIN
     >     ,PCO_MAX,PCO_LIM)
         !.. Check to adjust DTMAX and TWIDT near the end of run
         CALL ADJDTMAX(ITAU,DT,TSTOP,TWIDT,DTMAX)
        ENDIF

        !... set up the magnetic field line spatial grid .....
        IF(JTI.EQ.1.OR.L_STAG.NE.0) CALL RFIELD(7,PCO,RE,Z0,SCAL,SCALK,
     >    ZLBDY,Q,VPERPN,VTOT,COSDIP,JB,JR,JNQ,JTI,ZHI,ZLO,JVERT,JVERC)


        !... see if auroral file exists and adjust the vertical grid spacing
        CALL TFILE(14,AUR_FILE)
        IF(AUR_FILE.NE.0.AND.GRD.LT.0.5) GRD=1.0
        IF(AUR_FILE.NE.0.AND.Z0.GT.80.0) Z0=80.0

        !--- vertical altitude grid and initial densities for odd N and N2*
        CALL VERGRD(IVLPS,GRD)

        !--- evaluate the geographic location and magnetic declination for headers
        CALL AMBS(JVERT,Z(JVERT),GL(JVERT),D,T,CHIN,SATN,GLONN,GLATN)
        CALL AMBS(JVERC,Z(JVERC),GL(JVERC),D,T,CHIF,SATF,GLONF,GLATF)

        !....... write headers in the files on first time through
        IF(JTI.EQ.1) THEN
          CALL PHEAD(3,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
          CALL PHEAD(4,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
          CALL PHEAD(6,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
          CALL PHEAD(9,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
          CALL PHEAD(12,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
          CALL PHEAD(17,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
          CALL PHEAD(18,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
          CALL PHEAD(56,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
          CALL PHEAD(66,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))

          !-- Print the run control parameters - PCOFM obsolete
          CALL RUNPRN(6,ISKP,Z0,HPEQ,FPAS,HFAC,UVFAC,AUR_FILE,PCO)

          !--  Print info about data in wind file and ExB
          CALL PWINDS(IWIND,SATN,SATF,SEC,TSTOP,COSDIP(JK),ZTOPTE,
     >      NCOL_EXB,IODDN,IVIBN2)
          WRITE(6,116) ISOL,ITF,JMAX,IWIND,MAXIT,IVLPS,JTIMAX
     >    ,NINT(HPEQ),FPAS,HFAC,NINT(DTSET),NINT(DTMAX),NINT(TWIDT)
     >   ,TSTOP,HEPRAT,NINT(Z0),COLFAC,(ETRAN-43200.0)/86400.0
        ENDIF

 116    FORMAT(/1X,'ISOL ITF JMAX IWIND MAXIT IVLPS  JTIMAX '
     >   ,3X,'HPEQ',2X,'FPAS',1X,'HFAC',2X,'DTSET',2X,'DTMAX',2X
     >   ,'TWIDT',2X,'TSTOP',3X,'HEPRAT',5X,'Z0',2X,'COLFAC Eph_tran'
     >   ,/2I4,4I6,I9,I8,2F5.1,3I7,2F9.2,I5,F7.1,F7.2)

        CALL AMBS(JVERT,Z(JVERT),GL(JVERT),D,T,CHIN,SATN,GLONN,GLATN)
        CALL AMBS(JVERC,Z(JVERC),GL(JVERC),D,T,CHIF,SATF,GLONF,GLATF)

        ! When fitting topside Te, no plasmaspheric heating by e*
        IF(ZTOPTE.GT.0.0) FPAS=2.0

        IF(ITF.LE.-9) THEN
         WRITE(6,55) SEC/3600.,JDAY,(AP(IAP),IAP=1,7),NINT(F107)
     >      ,NINT(F107A)
         GO TO 221
        ENDIF

        !..  index of print lower bound for electron heating and cooling
        DO J=1,JMAX/2
           IF(Z(J).LT.ZLO) JZLO=J  
        ENDDO
     
 55     FORMAT(' UT=',F7.2,', DAY=',I8,'  AP(1..7)= ',7(1X,F5.1)
     >    ,' F10.7=',I4,' F10.7A=',I4)
      !.. Evaluate the magnetic field line cmpt. of neutral wind UN(J)
      !.. Alternative winds are read from file. Otherwise Hedin's winds
      !.... RATIN and RATIS are used to normalize NmF2 to IRI if IWIND<-1
      IF(IWIND.NE.0)
     >    CALL ALTWIN(DT,Z,N,ISOL,IWIND,TIMSAV,TSTOP,IWTIM,IWW,JMAX
     >    ,COSDIP,UNN,UNC,HMF2N,HMF2S,NMF2N,NMF2S,UN,DNMF2N,DNMF2S
     >    ,DHMF2N,DHMF2S,RATIN,RATIS,TOPTEN,TOPTES,ZTOPTE,APFILE)
     
        !---------- initialization and neutral atmosphere loop ----------------
        DO J=1,JMAX
           !..... o+ + n2* reaction rate factor .....
           VC(J)=1.0
           IF(RCON(J).LT.1.0) RCON(J)=1.0
           IF(RCON(J).GE.1.0) VC(J)=RCON(J)
           IF(JTI.EQ.1.AND.GL(J)*57.2958.GE.BLATPR) JBLATPR=J  !.. 2020-08-30
           IF(J.LE.JQ.AND.Z(J).LT.250.0) J250=J     ! index for F2 region

           !... Obtain neutral densities -- AMBS calls MSIS ......
           CALL AMBS(J,Z(J),GL(J),D,T,CHI,SAT,GLON,GLATJ)
           HE(J)=D(1)
           ON(J)=D(2)
           HN(J)=D(7)
           N2N(J)=D(3)
           O2N(J)=D(4)
           TN(J)=T(2)
           SZA(J)=CHI
           N4S(J)=D(8)
           !... Get GW data to modify winds. 2023-02-06
           !... GWDENTN = TN,OX,O2,N2,HE,H perturbations from GWs   
           !.. Gravity wave horiz component GW_BMER added 2023-02-06
           CALL GRAVITY_WAVE_DATA
     >         (0,J,SEC/3600.0,REAL(Z(J)),GLATN,GWDENTN,GWIND,GW_BMER)

           !.. If wind is not from HMW14, add in Gravity wave wind
           IF(IWIND.NE.0) THEN
              UN(J)= UN(J) -ABS(COSDIP(J))* GW_BMER * 100.0
              IF(J.EQ.JK) UNN= UNN + GW_BMER  *100.0
              IF(J.EQ.JC) UNC= UNC - GW_BMER  *100.0
              BCOMPU=0.0  !.. Zero for the write statement below
           ENDIF

           !.. from Hedin's model (BCOMPU=component along magnetic meridian)
           IF(IWIND.EQ.0) THEN
              CALL HEDWIN(J,JDAY,SX,REAL(Z(J)),GLATJ,GLON,SAT,F107A
     >          ,F107,AP,BCOMPU,UHED)
              UN(J)= -ABS(COSDIP(J))* (BCOMPU+GW_BMER) * 100.0
              IF(J.EQ.JK) UNN= (BCOMPU+GW_BMER) * 100.0
             IF(J.EQ.JC) UNC=-(BCOMPU+GW_BMER) * 100.0
           ENDIF

           !.. Output added 2023-02-09
           IF(J.EQ.JK) THEN
             IF(ABS(GWDENTN(2))+ABS(GWDENTN(3))+ABS(GWBMER).GT.0.0)
     >       WRITE(6,'(A,A,/2X,8F8.2,1P,22E10.2)') 
     >      '      UT    ALT   HWM_Bmer GW_Bmer   GWNth   GWEst   GWVer'
     >      ,'   GW_TN   GW_OX     GW_O2     GW_N2      GW_H     GW_HE'
     >      ,SEC/3600.0,Z(J),BCOMPU,GW_BMER,GWIND(1),GWIND(2),GWIND(3)
     >      ,(GWDENTN(I),I=1,6)
            ENDIF     
            !..... Calculate centrifugal force and add to gravity
            !..      IF(JTI.EQ.1) CALL CENFOR(IEXP,BLON,GL(J)*57.296,Z(J),CF)
            !..      IF(JTI.EQ.1.AND.IVLPS.GE.0) GR(2,J)=GR(2,J)+CF

            !..... to set up initial temperature and density profiles
            IF(HPEQ.GT.0.0) CALL PROFIN(ISOL,J,JMAX,Z(J),F107,N,TI,TN
     >       ,SL,GR,ON,HN,HPEQ,HEPRAT,XIONN,IVLPS)

            !..... printing of ambient parameters.................
            IF(NINT(SEC/3600.).NE.-43) GO TO 17
              JNEUT=JNEUT+1
              IF(JNEUT.EQ.1) WRITE(3,112)
 112        FORMAT(/'  J      ALT    GL     SZA     GR    TN     N(O)'
     >       ,5X,'N(H)     N(N2)     N(O2)     N(HE)   WIND')
            IF(J.LT.JMAX/2.AND.Z(J).GE.100.AND.Z(J).LE.350)
     >         WRITE(3,'(I4,F8.0,2F7.2,2F7.0,1P,5E9.1,0P,F9.1)')
     >       J,Z(J),GL(J),SZA(J),GR(2,J),TN(J),ON(J),HN(J)
     >      ,N2N(J),O2N(J),HE(J),UN(J)/100.0
    
 17        CONTINUE
      ENDDO
      
C ---------------- end neutral atmosphere loop -------------------

        IF(HPEQ.GT.0) HPEQ=-HPEQ
        !.. Obtain symmetric ambient conditions - diagnostics
        IF(ISYM.NE.0) CALL SYM(JMAX,HE,EHT,UN,Z,SZA)

        !--- transfer vertical grid odd N to field line grid. If no odd N,
        !---- a flag is set in TRSTR 
        CALL TRSTR(IVLPS)

        !.... evaluate primary production rates from primpr
        DO J=1,JMAX
          IF(J.LE.JMAX/2) EFAC=EFACN
          IF(J.GT.JMAX/2) EFAC=EFACS
          CALL PRIMPR(J,Z,ON,N2N,O2N,HE,SZA,TN,JMAX,SETCHI,IVLPS,EFAC)
        ENDDO
        !..  pe2s=2-stream prog to get electron heating rate;
        DO I=1,3
        DO J=1,JMAX
          EHT(I,J)=0.0
        ENDDO
        ENDDO
        IF(CHIN.LE.2.0.OR.CHIF.LE.2.0) CALL PE2S(Z,ON,N2N,O2N,HE,BM
     >   ,BG,JMAX,DH,N,TI,EHT,TN,SZA,FPAS,IVLPS,-1.0E22,ISOL,SL)

        !..... If there is an EUV flux file need to reevaluate primary rates
        IF(FEXISTS.EQ.-9) CALL TFILE(42,FEXISTS)
        IF(FEXISTS.NE.0) THEN
          DO J=1,JMAX
            IF(J.LE.JMAX/2) EFAC=EFACN
            IF(J.GT.JMAX/2) EFAC=EFACS
            CALL PRIMPR(J,Z,ON,N2N,O2N,HE,SZA,TN,JMAX,SETCHI,IVLPS,EFAC)
          ENDDO
        ENDIF

        !.... AURSUB calculates auroral electron production rates
        UTHRS=SEC/3600.
        CALL AURSUB(1,JMAX/2,UTHRS,Z,ON,O2N,N2N,N,TI,EHT,DTMAX,AURALT)

        !-- Sum the EUV, photoelectron, and auroral production rates
        CALL SUMPRD(JMAX)

        !.... adjust the plasmaspheric heating rate
        CALL CAL_EHEAT(ISOL,JQ,JHN,JHC,JHX,JMAX,JK,JC,ZTOPTE
     >    ,Z,TN,TI,PCO,HFAC,TIMEON,TIMEOF,F107,TOPTEN,TOPTES
     >    ,SEC,SATN,SATF,SATEQ,TISAV,N,EHT,HEATF)

        !.. loop to calculate ionization rates
        DO J=1,JMAX 
          PHION(J)=1.0E-22
          N(3,J)=1.0E-22
          DO IX=1,9
            FD(IX)=0.0
          ENDDO
          IF(Z(J).GE.Z0.AND.Z(J).LE.ZPROD) THEN
            !.. CALL cminor to get no+, o+(2d), n2+ & o+(2p) densities,
            CALL CMINOR(0,J,0,IPR,FD,7,HE(J))
            !-- PHION=total O+(4S) prod, including EUV, e*, dissoc of O2 and
            !-- minor ions. BCKPRD = small background production to eliminate 
            !-- instability below 200 km
            BCKPRD=2.0E-10*N2N(J)*EXP(-5.0E-14*N2N(J)*TN(J))
            PHION(J)=SUMION(1,7,J)+SUMION(2,4,J)+SUMION(2,5,J)+FD(9)
            SUMION(3,9,J)=FD(9)   !%%%%%%%%%
            PHION(J)=PHION(J)+BCKPRD
          ENDIF
          !.. Sum minor ions N+, NO+, O2+ for electron density at low altitudes
          !.. Leave out N2+, O+(2P), O+(2D) causes [e] gradient problems
          N(3,J)=XIONN(1,J)+FD(1)+FD(2)+(FD(3)+FD(4)+FD(5))*0.0000
        ENDDO

        !... Obtain symmetric ambient conditions - diagnostics
        IF(ISYM.NE.0) CALL SYM(JMAX,HE,EHT,UN,Z,SZA)

        !----- O+, H+, DENSITIES + electron and ion TEMPERATURES ---
        IF(I1.GT.0)  CALL LOOPS(S,NEQ,N,TI,ITAU,JB,SEC,MAXIT,ISOL
     >    ,JK,JC,IWIND,RATIN,RATIS,HE)

        !... testing for nonconvergence
        IF(ITER.EQ.-1) CALL RUN_ERROR    !.. print ERROR warning in output files
        IF(ITER.EQ.-1) STOP 1        !.. could not find a solution

        IF(I1.GT.0.AND.ITER.LE.-2) GO TO 221  !.. jump to end of time loop

        !.. Obtain symmetric ambient conditions - diagnostics
        IF(ISYM.NE.0) CALL SYM(JMAX,HE,EHT,UN,Z,SZA)

        !....... diagnostic print  from subroutines TFIJ and DFIJ ....
        !-- Moved here 21 April 1992 to ensure continuity equation balance
        IF((REAL(ITAU)/3600.0.GE.TSTOP-0.5).AND.JWR.GT.0) THEN
          WRITE(6,*) ' Writing diagnostic info in FOR03*.DAT files'
          IF(JWR.EQ.2) CALL AVDEN(JMAX1,TI,4,ZLO,ZHI)
          DO J=2,JMAX1
            IF(JWR.EQ.1.AND.Z(J).GT.ZLO.AND.Z(J).LT.ZHI)
     >        CALL TFIJ(J,4,IPR,N,TI,F,JB,JBS,V,HE)
            IF(JWR.EQ.2.AND.Z(J).GT.ZLO.AND.Z(J).LT.ZHI)
     >        CALL DFIJ(J,4,IPR,N,TI,F,JB,JBS,V)
          ENDDO
        ENDIF  !.. end print

        !--------------------------- He+ DENSITIES --------------------
c       IF(IABS(IVLPS).GE.9)
        CALL XION(S,NEQ,N,TI,ITAU,JB,SEC,MAXIT,ISOL,9)

        !---------------------------  N+ DENSITIES --------------------
c       IF(IABS(IVLPS).GE.11)
        CALL XION(S,NEQ,N,TI,ITAU,JB,SEC,MAXIT,ISOL,11)

        !... testing for nonconvergence
        IF(ITER.EQ.-1) CALL RUN_ERROR  !.. print ERROR warning in output files
        IF(ITER.EQ.-1) STOP 1        !.. could not find a solution
        IF(I1.GT.0.AND.ITER.LE.-2) GO TO 221

        !---------- N(2D), N(4S), NO, and N2* DENSITIES -----------
        !----- The odd N and N2* time step is never less than TWIDT
        VDT=SEC-SECSAV
        IF(ISOL.EQ.-1.AND.VDT.LT.TWIDT.AND.ITF.EQ.1) GO TO 277
        IF(VDT.LE.0) VDT=DT
        SECSAV=SEC

        IF(IODDN.GE.1.OR.IVIBN2.GE.1)
     >    CALL VERCON(J250,NOCON,S,JMAX,999,VDT,GL(JVERT),GRD,ISOL,IVLPS
     >   ,IVIBN2,HMF2N,HMF2S,SATN,SATF,CHIN,CHIF)

        !... Resetting densities to previous values due to odd N or N2*
        !... non-convergence
        IF(NOCON.EQ.0)  THEN
           CALL SAVDEN(1)       !... save good densities
        ELSE
           CALL RECDEN(ITER)    !... recall previous densities
           GO TO 221
        ENDIF

        IF(ISOL.EQ.0) ITF=ITFS
        ISOL=-1

 277    CONTINUE

        ! ******************** STORE RESULTS *********************
        !....... write a new file at each time step ........
        CALL WRITE_BQ(1,JMAX,SCALK,Q)
        CALL WRSTRT(PCO,Z0,SCAL,GRD,GLONN,GLONF,GLATN,GLATF,UNN,UNC
     >   ,HMF2N,HMF2S,NMF2N,NMF2S,PLAT,PLON)

        !-- Printing summary densities and temperatures in ASCII files
        !-- First, write warning about initial conditions then data in file
        IF(ICAUTN.EQ.0.AND.REAL(ITAU)/3600.0.GT.12) CALL PRCAUT(ICAUTN)
        !.. Get the winds at altitude ZLO in both hemispheres
        CALL HEDWIN(JZLO,JDAY,SX,REAL(ZLO),GLATN,GLONN,SAT,F107A,F107,AP,
     >    BCOMPU,UHED)                 !.. North wind
        UMZ(1,1)=UHED(1)
        UMZ(1,2)=UHED(2)
        JZLOC=JMAX+1-JZLO
        CALL HEDWIN(JZLOC,JDAY,SX,REAL(ZLO),GLATF,GLONF,SAT,F107A,F107,
     >    AP,BCOMPU,UHED)                 !.. South wind
        UMZ(2,1)=UHED(1)
        UMZ(2,2)=UHED(2)
        CALL DIAGPR(ISOL,JDAY,JB,JK,JC,JQ,JNQ,JFQ,HMF2N,HMF2S,UMZ,NMF2N
     >  ,NMF2S,SATN,SATF,CHIN,CHIF,UNN,UNC,VTOT,PCO,VPERPN,COSDIP)

        !-- Output for comparison with IRI model 20 APR 1992
        CALL EDEN170(JMAX,Z,N,J170,N170KMN,N170KMS)  !.. get [e] at 170 km
        CALL IRIPRN(JVERT,JVERC,IRIHED,SEC,SATN,CHIN,NMF2N
     >    ,DNMF2N,HMF2N,DHMF2N,SATF,CHIF,NMF2S,DNMF2S,HMF2S,DHMF2S
     >    ,IDAY,GLATN,GLONN,GLATF,GLONF,F107A,N170KMN,N170KMS)

        !... output transition altitudes 24/2/2000
        CALL TRAN_ALT(N,JMAX,Z,ON,HN,SATN,SATF,UTHRS)

        !-- Output production frequencies
        DO J=1,JMAX/2
          IF(Z(J).LT.600) JFREQN=J
        ENDDO
        JFREQS=JMAX+1-JFREQN
        CALL EUVPRN(JK,JC,JFREQN,JFREQS,SEC,SATN,CHIN,SATF,CHIF,IDAY
     >     ,F107,F107A,TI,Z)
        CALL CAL_STATS(SEC,DT,HMF2N,DHMF2N,HMF2S,DHMF2S)

        !.. Calculate Richards and Torr limiting H+ flux 
        CALL FLXLIMIT(JR,JQ,JK,Z,N,TI,BM,GR,ON,HN,TN,SL,JZ,FYLIM,HPLOSS,
     >    SH_OPLUS,ZTRAN)
        FLIPFLX=N(2,JNQ)*U(2,JNQ)*BM(1)/BM(JNQ)

         !.. Find the peak H+ density to use as transition altitude
        JHP=JQ   !.. begin at hmF2 and go up the field line
        DO J=JQ-1,JK,-1
          !.. Check for peak density altitude
          IF(N(2,J).GT.N(2,J-1).AND.N(2,J).GT.N(2,J+1))
     >      JHP=J
        ENDDO
      
         IF(MOD(JTI-1,18).EQ.0) WRITE(6,317)
 317     FORMAT(/10X,'[------Northern hemisphere-------]',5X,
     >   '[------Southern hemisphere------]',/3X,'  UT    LT  SZA'
     >   ,3X,'NmF2   HmF2  N2* Lat Lon     LT  SZA   NmF2   HmF2'
     >   ,2X,'N2* Lat Lon  Kp   F107   JTI  DelT ITR   L-Sh   YYYYddd')
 117     FORMAT(F8.2,F6.2,I4,1P,E9.2,I5,0P,F5.2,2I4,2X,F6.2,I4,1P
     >   ,E9.2,I5,0P,F5.2,2I4,F5.1,I5,I8,I5,I4,F9.3,I8)
         WRITE(6,117) (SEC)/3600.,SATN,NINT(57.3*CHIN),NMF2N
     >   ,NINT(HMF2N),RCON(JK),NINT(GLATN),NINT(GLONN),SATF
     >   ,NINT(57.3*CHIF),NMF2S,NINT(HMF2S),RCON(JC),NINT(GLATF)
     >   ,NINT(GLONF),KPREAL,NINT(F107),JTI,NINT(DT),ITER,PCO,JDAY

        !.. Special output above an equatorial station 2020-08-30
        !.. Only output if BLATPR is not exactly zero
        IF(JBLATPR.GE.1.AND.JBLATPR.LE.JMAX.AND.ABS(BLATPR).GT.0.001) THEN
          IF(JTI.EQ.1)
     >      CALL PHEAD(43,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
          CALL AMBS(JBLATPR,Z(JBLATPR),GL(JBLATPR),D,T,CHIX,SATX,
     >       GLONX,GLATX)
          IF(JTI.EQ.1) WRITE(43,'(/5X,A,A,2F5.1,A,F5.1,F6.1,A)') 
     >      ' Ionosphere quantities located on field line above',
     >      ' (MagLat. MagLon.)=(', BLATPR,BLON,'); [GeoLat. Geolon]=['
     >       ,GLATX,GLONX,']'
          IF(JTI.EQ.1) WRITE(43,'(11X,A,/3X,A,7X,A)') 
     >     'Ne = electron density; min+ = O2+ + NO+ + N2+ + N+'
     >  ,'UT       LT      SZA      ALT     Tn      Ti      Te       Ne'
     >  ,'O+       H+      min+     He+       OX      O2       N2'

         WRITE(43,'(3F8.2,F10.2,3F8.2,1P,9E9.2)') 
     >     SEC/3600.,SATX,57.3*CHIX
     >    ,Z(JBLATPR),TN(JBLATPR),TI(1,JBLATPR),TI(3,JBLATPR)
     >    ,N(1,JBLATPR)+N(2,JBLATPR)+N(3,JBLATPR)+N(4,JBLATPR)
     >    ,N(1,JBLATPR),N(2,JBLATPR),N(3,JBLATPR),N(4,JBLATPR)
     >    ,ON(JBLATPR),O2N(JBLATPR),N2N(JBLATPR)
        ENDIF !.. End special output 2020-08-30
      
        !.... Write topside heat flow and Brace and Theis Te and write
        IF(JHFLX.EQ.0) THEN
          !..JJ=JK    !.. altitude for cooling rates. Can use JC also
          JJ=JZLO  !.. altitude for cooling rates. Can use JK or JC also
          WRITE(56,893) NINT(Z(JHX)),NINT(Z(JJ))
 893      FORMAT(/3X,'Hfac= heating factor, heat_eq= heating rate at'
     >    ,1X,'equator, FLUX= heat flux at',I5,' km (eV cm-2 s-1)'
     >    ,3X,'other heating and cooling rates are at',I7,'km'/
     >    ,3X,'The last 5 colums are for the Richards and Torr 1985'
     >    ,1X,'North limiting H+ flux, HPLOSS = H+ loss above Z0'/)
          IF(NINT(ZTOPTE).GT.200) THEN
            WRITE(56,894) NINT(Z(JHN)),NINT(Z(JHN))
 894        FORMAT(' * Te data constraint imposed at ',I5,' km'
     >           ,1X,' and Tcalc is also at',I5,' km')
          ELSE IF(NINT(ZTOPTE).EQ.1) THEN
             WRITE(56,895)
 895        FORMAT(' * Te data constraint imposed at hmF2'
     >           ,1X,' and Tcalc is also at hmF2')
          ELSE
            WRITE(56,896) NINT(Z(JHX))
 896        FORMAT(' * NO Te data constraint, Tcalc is at',I5,' km')
          ENDIF
        ENDIF    !... end header info


        IF(SATN-SATNSAV.LT.0)  WRITE(56,897)
 897      FORMAT(2X,'  UT      LT     hmF2    SZA Tmeas  Tcalc Ti_cal'
     >    ,2X,'Hfac    FLUX    heat_eq  Tn_hm Ti_hm   e_heat'
     >    ,2X,'ion_cool  N2_cool  O_cool   e_N2D      Ne     H_C'
     >    ,11X,'Z0       ZTRAN   PhiLIM    PhiFLIP   HPLOSS    PhiRAT'
     >    ,5X,'Oplus    Hy_Z0    SH_OPLUS ZHpeak   Tdiff')
        SATNSAV=SATN

        IF(NINT(ZTOPTE).LT.1) JHN=JHX     !... no Te constraint
        JHFLX=1
        HFLUX1=7.7E+5*TI(3,JHX)**2.5*(TI(3,JHX)-TI(3,JHX-1))/
     >      (SL(JHX)-SL(JHX-1))
        !...  DNE=N(1,JHN)+N(2,JHN)+N(3,JHN)+N(4,JHN)
        !... CALL BRACE(Z(JHN),DNE,TI(3,JHN),TB,TEOTB)

        CALL ECRATS(0,JJ,Z(JJ),TI(3,JJ),TI(2,JJ),TN(JJ),N(1,JJ)
     >     ,N(2,JJ),N(3,JJ),ON(JJ),N2N(JJ),O2N(JJ),HN(JJ),HE(JJ)
     >     ,EHT(3,JJ),TLSS,FOL,COOL)
           CALL CMINOR(0,JJ,0,IPR,FD,7,HE(JJ))

        WRITE(56,218) SEC/3600,SATN,Z(JK),NINT(57.3*CHIN)
     >      ,NINT(TOPTEN),NINT(TI(3,JHN)),NINT(TI(1,JHN)),HEATF
     >      ,HFLUX1,EHT(3,JQ),NINT(TN(JJ)),NINT(TI(1,JJ))
     >      ,(COOL(K),K=1,4),FD(8)*2.4,N(1,JJ)+N(2,JJ)+N(3,JJ)
     >      ,COOL(1)+FD(8)*2.4-COOL(5),Z(JZ),ZTRAN,FYLIM
     >      ,FLIPFLX,HPLOSS,FLIPFLX/FYLIM,N(1,JZ),HN(JZ),SH_OPLUS/1E5
     >      ,Z(JHP),(TOPTEN-TI(3,JHN))

 218   FORMAT(3F8.2,4I6,F8.2,1P,E10.2,1P,E10.2,2I6,1P,E10.2,5E9.2
     >   ,E10.2,0P,2F10.1,1P,7E10.2,0P,2F8.0)

        !... for printing ion-neutral cooling. not ready yet
        !..        CALL ION_NEUTRAL(0,N(1,JK),N(2,JK),ON(JK),N2N(JK),O2N(JK)
        !..     >   ,HN(JK),HE(JK),TN(JK),TI(2,JK),Z(JK),COOL)

        !.... set a switch to prevent writing into unit 8 until midnight if
        !.... setting up new field line
        IF(ISOL.EQ.0) WUNIT8=0
        !...      IF(ISOL.NE.0.AND.SEC.GE.0.0) WUNIT8=1

        !--- CALL subroutine to write main data on unit 8 at intervals TWIDT
        !--- Data is not written until run has settled down
        IF(ITF.NE.1) GO TO 221
        IF(WUNIT8.NE.0.OR.SATN.LT.2.OR.SATF.LE.2)  THEN
          IF(SEC-SECWF8.GE.TWIDT.AND.NOWR8.EQ.0) THEN
            CALL WF8DAT(IVLPS,JQ,PCO,Z0,SCAL,GRD,GLONN,GLONF,GLATN,
     >       GLATF,ZLBDY,FPAS,UV,HFAC,TIMEON,TIMEOF,AURALT,PLAT,PLON)
            CALL WF8D2(JK,JC,UNN,UNC,IVLPS,JQ,UE)
            SECWF8=SEC
          ENDIF
          WUNIT8=1  !... switches on print from here on
        ENDIF
        ! ******************** END STORE RESULTS *************************

        !... print densities and temperatures from subroutine loops....
        IF(JTIMAX.NE.JTI) GO TO 221
        IF(ITF.EQ.3.OR.ITER.LE.0.OR.JTIMAX.EQ.1) ITSAV=-1000000
        IF(ITAU-ITSAV.LT.3600) GO TO 221
        ITSAV=ITAU
        IF(IPP.EQ.0.OR.JWR.EQ.0) GO TO 221
        WRITE(3,859)
        WRITE(9,859)
 859    FORMAT(/'  J ',3X,'ALT',3X,'TE',4X,'TI',3X,'(O+)',5X,'(H+)',6X
     >   ,'VO+',6X,'VH+',5X,'(He+)',5X,'VHe+',5X,'EHT    (N+)     VN+')
        DO J=1,JQ
          JU=JMAX+1-J
          IF(Z(J).LT.ZLO.OR.Z(J).GT.ZHI) GO TO 854
          IF(J.EQ.1.OR.J.EQ.JQ) GO TO 852
          IF(MOD(J,IPP).EQ.0) GO TO 852
          GO TO 854
 852      CONTINUE

          NE=N(1,J)+N(2,J)+N(3,J)
          CALL BRACE(Z(J),NE,TI(3,J),TB,TEOTB)
          WRITE(3,888) J,NINT(Z(J)),NINT(TI(3,J)),NINT(TI(1,J))
     >     ,N(1,J),N(2,J),U(1,J),U(2,J),N(4,J),UE(1,J),EHT(3,J)
     >     ,XIONN(1,J),XIONV(1,J)
          NE=N(1,JU)+N(2,JU)+N(3,JU)

          CALL BRACE(Z(JU),NE,TI(3,JU),TB,TEOTB)
          WRITE(9,888) JU,NINT(Z(JU)),NINT(TI(3,JU)),NINT(TI(1,JU))
     >     ,N(1,JU),N(2,JU),U(1,JU),U(2,JU),N(4,JU),UE(1,JU),EHT(3,JU)
     >     ,XIONN(1,JU),XIONV(1,JU)
 854      CONTINUE
        ENDDO
 221    CONTINUE
      ENDDO   !.. End of main loop
      
 64   FORMAT(2I6,3D12.4,/11F10.2)
 65   FORMAT(16A8)
 66   FORMAT(32A4)
 68   FORMAT(22I6)
 111  FORMAT(I4,F8.0,2F7.2,2F7.0,1P,22E9.1)
 167  FORMAT(10E10.2)
 888  FORMAT(I4,I6,I6,I5,1P,2E9.2,1P,5E9.1,1P,12E9.1)
      END

C::::::::::::::::::::::: CAL_EHEAT ::::::::::::::::::::::::::::
      !.. Adjusting plasmaspheric heating rate.
      SUBROUTINE CAL_EHEAT(ISOL,JQ,JHN,JHC,JHX,JMAX,JK,JC,ZTOPTE
     >   ,Z,TN,TI,PCO,HFAC,TIMEON,TIMEOF,F107,TOPTEN,TOPTES,SEC
     >   ,SATN,SATS,SATEQ,TISAV,N,EHT,HEATF)
      IMPLICIT NONE
      INTEGER ISOL,JQ,JHN,JHC,JHX,JMAX,IL,JC,JK,J,IH
      REAL ZTOPTE,ZHFLX,TIMEON,TIMEOF,F107,SATN,SATS,PROP,SATEQ
     >   ,TMAX,SIGT,LMAX,SIGL,FAC_RING,RING_CURR,AMP_RING
     >   ,RC_TIME_CON,SEC
      DOUBLE PRECISION Z(401),TN(401),TI(3,401),HFAC,EHT(3,401),PCO,
     >  TOPTEN,TOPTES,HEATF,TISAV(3,401),N(4,401),HEATFS,TOPTENSV,TETEMP
      !... parameters for the ring current model
      COMMON/RC/AMP_RING,TMAX,SIGT,LMAX,SIGL,RC_TIME_CON,PROP
      DATA ZHFLX/486.0/     !.. altitude of topside heat flux
      DATA HEATFS/0.0/
      
      !.. First determine grid point where Te normalization will be done
      JHN=JMAX/2
      JHC=JMAX/2
      JHX=JMAX/2
      DO 545 IL=2,JMAX/2
        IF(ZTOPTE.LT.(Z(IL+1)+Z(IL))/2.AND.ZTOPTE.GE.(Z(IL)+Z(IL-1))/2)
     >     JHN=IL
        IF(ZHFLX.LT.(Z(IL+1)+Z(IL))/2.AND.ZHFLX.GE.(Z(IL)+Z(IL-1))/2)
     >     JHX=IL
 545  CONTINUE

      JHC=JMAX+1-JHN            !..... Conjugate index
      !... reset if topside Te is at hmF2
      IF(JK.GT.1.AND.NINT(ZTOPTE).EQ.1) JHN=JK
      IF(JK.GT.1.AND.NINT(ZTOPTE).EQ.1) JHC=JC

      !... determine factor for ring current heating of ions
      IF(AMP_RING.GT.0.0.AND.RC_TIME_CON.GT.0.0) THEN
        ZTOPTE=-1
        FAC_RING=AMP_RING*RING_CURR(SATEQ,REAL(PCO),TMAX,SIGT,LMAX,SIGL)
        FAC_RING=FAC_RING*EXP(-SEC/3600.0/RC_TIME_CON)   !.. time decay
        !..WRITE(6,95) NINT(TI(1,JQ)),NINT(TI(3,JQ))
        !.. >     ,NINT(TN(JQ)),FAC_RING,FAC_RING*(N(2,JQ)+N(1,JQ))
c 95     FORMAT(' RC ',3I9,1P,2E11.3) 
        DO J=1,JMAX
          IF(Z(J).GT.2000) THEN
             EHT(1,J)=FAC_RING*(1.-PROP)*(N(2,J)+N(1,J))+EHT(1,J)
             EHT(3,J)=FAC_RING*PROP*(N(2,J)+N(1,J))+EHT(3,J)
          ENDIF
        ENDDO
      ENDIF
      
      IH=0  !.. switch to turn on heating
      !.. before midnight
      IF(HFAC.GT.0.01.AND.SATEQ.GE.TIMEON.AND.SATEQ.LE.TIMEOF) IH=1
      !.. after midnight assuming TIMEOF >= 24
      IF(HFAC.GT.0.01.AND.TIMEOF.GE.24.AND.SATEQ+24.LE.TIMEOF) IH=1
      IF(IH.EQ.1) THEN
        !... Unknown heat source. HFAC=1 -> HF=1.0E9 eV cm-3s-1 OR 
        !... if <0 a straight multiplicative factor at all altitudes
         DO J=1,JMAX
           IF(Z(J).GT.1000)  EHT(3,J)=2.5*HFAC/PCO**4 + EHT(3,J)
         ENDDO
         ZTOPTE=-1  !... signal to return
      ELSE
        !.... If HFAC < 0.0 use constant scaling
        IF(HFAC.LT.-0.01) THEN
          DO J=1,JMAX
           IF(ZTOPTE.GT.200.AND.Z(J).LT.ZTOPTE) EHT(3,J)=-HFAC*EHT(3,J)
           IF(ZTOPTE.LE.200) EHT(3,J)=-HFAC*EHT(3,J)
          ENDDO
           ZTOPTE=-1  !... signal to return
        ENDIF
      ENDIF

      IF(NINT(ZTOPTE).LE.0) RETURN
      IF(ISOL.NE.-1) RETURN  !.. only for regular run

      !.. Increased plasmaspheric heating. HEATF is used when reproducing
      !.. Te input from a file using heat flux algorithm
      !.. 2020-09-27: Added HEATFS and TOPTENSV. Also since the algorithm 
      !.. needs a Te difference, a small percentage is added to TOPTEN
      IF(ZTOPTE.GT.0) THEN
         !... HEATFS, TETEMP, and TOPTENSV smooth out oscillations in HEATF
         HEATFS=HEATF
         IF(TOPTEN.GT.TN(JHN)) 
     >     TETEMP=0.333*TOPTENSV+0.666*TOPTEN+0.04*ABS(TOPTEN-TN(JHN))
         IF(TOPTEN.GT.TN(JK)) HEATF=5E-8*(TETEMP**2.5-TI(3,JHN)**2.5)
         IF(TOPTES.GT.TN(JC)) HEATF=5E-8*(TOPTES**2.5-TI(3,JHC)**2.5)

         !.. make sure heating factor is reasonable. It can be negative in 
         !.. case a previous HEATF was too big.
         IF(HEATF.GT.15.0) HEATF=15.0
         IF(HEATF.LT.-2.0) HEATF=-2.0  

         HEATF=(HEATF+HEATFS)/2.0      !.. Smooth heat flux - not needed
         !.. If MSIS Tn and TOPTEN don't match could get weird profiles 
         !.. and non-convergence
         IF(HEATF.LT.0.0.AND.TI(3,JQ).LT.TETEMP+100) HEATF=0

         TOPTENSV=TETEMP   !.. save for smoothing
         
         !... Now adjust electron EHT(3,J) and ion heating rates
         DO J=1,JMAX
           IF(Z(J).GT.ZTOPTE+500) THEN
              EHT(1,J)=9.0*(1.-PROP)*HEATF/PCO**4 + EHT(1,J)
              EHT(3,J)=9.0*PROP*HEATF/PCO**4 + EHT(3,J)
           ENDIF
         ENDDO
      ENDIF
      
      RETURN
      END
C::::::::::::::::::::::::::: RING_CURR ::::::::::
C... This function fits the ring current using two Guassian distributions
      REAL FUNCTION RING_CURR(LT,       !.. local time
     >                         L,       !.. L-value
     >                      TMAX,       !.. time of max ring current
     >                      SIGT,       !.. width of time Guassian
     >                      LMAX,       !.. L-value of max ring current
     >                      SIGL)       !.. width of L Guassian
      IMPLICIT NONE
      REAL LT,L,TMAX,SIGT,LMAX,SIGL,DELTIME
      !... First term is for time. Second term is for L
      DELTIME=ABS(LT-TMAX)
      IF(ABS(LT-TMAX).GT.ABS(LT+24-TMAX)) DELTIME=ABS(LT+24-TMAX)
      RING_CURR=EXP(-(DELTIME/SIGT)**2) * EXP(-((L-LMAX)/SIGL)**2)
      RETURN
      END
C:::::::::::::::::::::: BLOCK DATA :::::::::::::::::::::::::::::
C..... This BLOCK DATA is used to initialize the arrays in the FLIP
C..... model. This is necessary for the IBM machine
      BLOCK DATA BLO3
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER ION
      REAL SEC,DEC,ETRAN,F107,F107A,AP,BLON,EUV,UVFAC,ZPAS,PASK,FPAS
      COMMON/PALT/ZLO,ZHI,JWR,JTI
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/LPS/SZA(401),EPSN,DC,ISPC,I1,I2,IPV,JTIMAX,IPP,ISKP
      COMMON/SOL/UVFAC(59),EUV
      COMMON/PANGS/ZPAS,PASK,IPAS,IPASC,FPAS
      DATA ZPAS /2000.0/, ZLO/0.0D0/,ZHI/0.0D0/,JWR/0/,JTI/0/
      DATA AP/7*0.0/,DEC/0.0/,ETRAN/0.0/,BLON/0.0/,F107/0.0/,F107A/0.0/
      DATA IDAY/0/,SEC/0.0/, DT/0.0D0/,DH/0.0D0/, ITF/0/,JMAX1/0/
      DATA IPV/0/,ITER/0/, JTIMAX/0/, EUV/0.0/, ISKP/1/, ISPC/5/
      DATA EPSN,  DC,I1,I2,IPP,    THF,EPS  ,ION, TF
     >   /  .999,.5D0, 1, 2, 4 , 1.0D0 ,1E-3, 2 ,1.0 /
      END

C::::::::::::::::::::::::::::::::: ACTUAL_DAY:::::::::::::::::::::::::
C...... Since IDAY does not change and UT continually increases, need to
C...... evaluate actual day (ID) and UT (SX) for MSIS, HWM, and IRI.
C...... Modified by P. Richards in November 2008 to allow for long runs
      SUBROUTINE ACTUAL_DAY(IDAY,    !.... FLIP start day YYYYDDD
     >                       SEC,    !.... Current UT in secs (continuous)
     >                        ID,    !.... actual day YYYYDDD (output) 
     >                        SX)    !.... UT (0-24) in secs (output)
      
      IMPLICIT NONE
      INTEGER SDIM                  !.. Dimension for array
      PARAMETER(SDIM=11)
      INTEGER I,IWR                 !.. Loop control, IWR= write control 
      INTEGER IDAY,ID               !.. FLIP start YYYYDDD, ID= current YYYYDDD
      INTEGER YSTART,DSTART         !.. Start year YYYY and start day DDD
      INTEGER NYRS,NDAYS,NSECS(SDIM)  !.. # of years, days, and seconds in these years
      INTEGER IDAYS                 !.. # Number of days in run so far
      REAL SEC,SX                   !.. UT in secs, actual UT in secs (output)
      DATA IWR/0/

      !.. Check the run is not too long
      IF(SEC.GT.SDIM*3.16E+7) THEN
         WRITE(6,55) SDIM
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF
 55   FORMAT(2X,'** ERROR: FLIP cannot do more than',I3,' year run **')

      YSTART=IDAY/1000               !.. Determine the start year YYYY
      DSTART=MOD(IDAY,1000)          !.. Determine the start day DDD

      !.. Determine the number of days and seconds simulated so far 
      !.. in order to determine the current year and the UT in that year
      NDAYS=0
      DO I=1,SDIM
        IF(MOD(YSTART+I-1,4).EQ.0) THEN     !.. leap year
           NDAYS=NDAYS+366
        ELSE
           NDAYS=NDAYS+365                  !.. normal year
        ENDIF

        !.. Time elapsed since the start day in seconds
        NSECS(I)=(NDAYS-DSTART+1)*24*3600

        !.. Jump out of loop when the time exceeds the current FLIP run time
        IF(NSECS(I).GT.INT(SEC)) THEN
           NYRS=I-1
           GOTO 10
        ENDIF
      ENDDO

 10   CONTINUE

      !.. Determine year and UT in the first year and then later years
      IF(NYRS.LE.0) THEN
        IDAYS=INT(SEC)/3600/24
        ID=(YSTART)*1000+DSTART+IDAYS
      ELSE                                   !.. after first year
        IDAYS=(INT(SEC)-NSECS(NYRS))/3600/24
        ID=(YSTART+NYRS)*1000+IDAYS+1
      ENDIF

      SX=AMOD(SEC,8.64E+4)  !.. UT 

      IF(MOD(ID/1000,4).EQ.0) THEN     !.. leap year
        IF(IDAYS.GT.366) THEN
          WRITE(6,*) '  ** ERROR > 366 days in leap year'
          CALL RUN_ERROR    !.. print ERROR warning in output files
        ENDIF
      ELSE
        IF(IDAYS.GT.365) THEN
          WRITE(6,*) '  ** ERROR > 365 days in normal year'
          CALL RUN_ERROR    !.. print ERROR warning in output files
        ENDIF
      ENDIF

      IWR=0  !..IWR+1  !.. IWR=0 switches off DEBUG write statements
      IF(IWR.EQ.1) WRITE(6,191) 
 191  FORMAT(5X,'IDAY',8X,'ID    NYRS IDAYS     SX',9X,'SEC',8X,'Days')
      IF(IWR.GT.0) WRITE(6,'(2I11,2I5,9I11)') IDAY,ID,NYRS,IDAYS,
     >    INT(SX),INT(SEC),NDAYS,NSECS(NYRS)

      RETURN
      END
C::::::::::::::::::::::::::::: TRAN_ALT :::::::::::::::::::::::::::::
C... This function calculates the transition altitudes from O+ to light
C... ions. Modified 2007-10-18 by P. Richards to do both hemispheres
      SUBROUTINE TRAN_ALT(N,    !.. ion densities
     >                 JMAX,    !.. number of grid points
     >                  ALT,    !.. altitude array
     >                   ON,    !.. oxygen array
     >                   HN,    !.. hydrogen array
     >                 SATN,    !.. Local Time in the Northern hemisphere
     >                 SATF,    !.. Local Time in the Southern hemisphere
     >                 UTHRS)   !.. Universal time
      IMPLICIT NONE
      INTEGER JTOPHPN,JTOPHPC,JTOPLIN,JTOPLIC  !.. transition indices
      INTEGER JMOLN,JMOLC                      !.. molecular transition indices
      INTEGER JMAX                        !.. number of grid points
      INTEGER J,JU,IWR                    !.. loop control variables
      REAL SATN,SATF,UTHRS                !.. Defined above
      REAL RAT1,RAT2                      !.. For calculating ratios
      DOUBLE PRECISION N(4,401),ALT(401)  !.. defined above
      DOUBLE PRECISION ON(401),HN(401)    !.. defined above
      
      !.. Determine the indices where the ion transitions take place
      DO 20 J=1,JMAX/2
        JU=JMAX+1-J                         !.. conjugate index
        IF(N(1,J).GT.N(2,J)) JTOPHPN=J
        IF(N(1,JU).GT.N(2,JU)) JTOPHPC=JU
        IF(N(1,J).GT.N(2,J)+N(4,J)) JTOPLIN=J
        IF(N(1,JU).GT.N(2,JU)+N(4,JU)) JTOPLIC=JU
        IF(ALT(J).LT.500.AND.N(1,J).LT.N(3,J)) JMOLN=J
        IF(ALT(JU).LT.500.AND.N(1,JU).LT.N(3,JU)) JMOLC=JU
 20   CONTINUE

      !------------ Northern Hemisphere --------------------
      !.. See if the next point will give a O+/H+ ratio closer to 1
      RAT1=N(1,JTOPHPN)/N(2,JTOPHPN) 
      RAT2=N(1,JTOPHPN+1)/N(2,JTOPHPN+1)
      IF(ABS(1-RAT1) > ABS(1-RAT2)) JTOPHPN=JTOPHPN+1

      !.. See if the next point will give a O+/(H+ + He+) ratio closer to 1
      RAT1=N(1,JTOPLIN)/(N(2,JTOPLIN)+N(4,JTOPLIN)) 
      RAT2=N(1,JTOPLIN+1)/(N(2,JTOPLIN+1)+N(4,JTOPLIN+1))
      IF(ABS(1-RAT1) > ABS(1-RAT2)) JTOPLIN=JTOPLIN+1

      !------------ Southern Hemisphere --------------------
      !.. See if the next point will give a O+/H+ ratio closer to 1
      RAT1=N(1,JTOPHPC)/N(2,JTOPHPC) 
      RAT2=N(1,JTOPHPC+1)/N(2,JTOPHPC+1)
      IF(ABS(1-RAT1) > ABS(1-RAT2)) JTOPHPC=JTOPHPC+1

      !.. See if the next point will give a O+/(H+ + He+) ratio closer to 1
      RAT1=N(1,JTOPLIC)/(N(2,JTOPLIC)+N(4,JTOPLIC)) 
      RAT2=N(1,JTOPLIC+1)/(N(2,JTOPLIC+1)+N(4,JTOPLIC+1))
      IF(ABS(1-RAT1) > ABS(1-RAT2)) JTOPLIC=JTOPLIC+1
     
      DATA IWR/0/  !.. IWR is used to write the header only once
      IWR=IWR+1
      IF(IWR.EQ.1) WRITE(18,95)
 95   FORMAT(/'   ---- approximate transition altitudes for O+ to H+'
     > ,1X,'and light ions (H+ + He+) and ion densities at those'
     > ,1X,'altitudes in both north and south hemispheres'
     > /,11X,'North Hemisphere O+ to H+ transition   North O+ to' 
     > ,1X,'(H+ + He+) transition',3X,'South Hemisphere O+ to H+'
     > ,1X,'transition    South O+ to (H+ + He+) transition'   
     > /,5X,'UT    LT     ALT      O+      H+      He+      ALT'
     > ,5X,'O+      H+       He+      LT     ALT     O+      H+'
     > ,7X,'He+      ALT     O+      H+       He+      OpMN   OpMS'
     > ,2X,'H/O_Nth  H/O_Sth') 

      WRITE(18,90) UTHRS,SATN,NINT(ALT(JTOPHPN)),N(1,JTOPHPN)
     > ,N(2,JTOPHPN),N(4,JTOPHPN),NINT(ALT(JTOPLIN)),N(1,JTOPLIN)
     > ,N(2,JTOPLIN),N(4,JTOPLIN)
     > ,SATF,NINT(ALT(JTOPHPC)),N(1,JTOPHPC),N(2,JTOPHPC),N(4,JTOPHPC)
     > ,NINT(ALT(JTOPLIC)),N(1,JTOPLIC),N(2,JTOPLIC),N(4,JTOPLIC),
     >  NINT(ALT(JMOLN)),NINT(ALT(JMOLC)),
     >  HN(JTOPHPN)/ON(JTOPHPN),HN(JTOPHPC)/ON(JTOPHPC)
 90   FORMAT(F8.2,F7.2,I7,1P,3E9.2,I7,1P,3E9.2,0P,F7.2,I7,1P,
     >    3E9.2,I7,1P,3E9.2,2I7,1P,2E9.1)
      RETURN
      END
C:::::::::::::::::: HP_KP ::::::::::::::::::::::::
C... Adjust the H+ density for high Kp.
C... Note that Binsack JGR 1967, page 5321 gives Lpp=6.0-0.6*Kp
C... Alternative from Moldwin et al. AGU Spring 2001 SM21A-11
C... Lpp(day)=4.7-0.31*Kp; Lpp(night)=5.5-0.37*Kp
C... Modified July 2008 to allow variable equatorial density instead of 300
      SUBROUTINE HP_KP(KP,      !... Kp index    
     >                PCO,      !... L-shell
     >              UTHRS,      !... Universal time
     >               JMAX,      !... max index on field line
     >                  Z,      !... altitude array
     >                  N)      !... ion densities
      INTEGER J,JMAX
      REAL UTHRS,ALPHA,KP,HP_DEP
      DOUBLE PRECISION PCO,N(4,401),Z(401)

      !.. equatorial density for a depleted flux tube. Note that plasma may be convected to 
      !.. larger L-shells but then some of it convected back from larger L-shells as activity abates
      !.. HP_DEP=30+1.22E4/PCO**4  !.. Old value replaced in August 2012 
      !.. HP_DEP=500.0-100.0*PCO   !.. gives a more gradual decreas with L
      !.. HP_DEP=240.0/PCO/PCO     !.. before 2022-12-20
      HP_DEP=3000.0/PCO/PCO        !.. Updated 2022-12-20
      IF(HP_DEP.LT.5.0) HP_DEP=5.0   !.. set floor value Chappell 1972, fig.4


      !.. No need to deplete already depleted flux tube
      IF(N(2,(JMAX+1)/2).LT.HP_DEP) RETURN

      ALPHA=HP_DEP/N(2,(JMAX+1)/2)   !.. reduction fac. at equator

      !.. If Kp is high enough, reduce H+ and He+ densities. Note that 24 hrs
      !.. must have elapsed since previous reduction, otherwise too unstable. 
      IF(KP.GE.((6.0-PCO)/0.6).AND.UTHRS-THPEQ.GT.24.) THEN !..Binsack
        DO J=1,JMAX
          N(2,J)=N(2,J)*ALPHA
          N(4,J)=N(4,J)*ALPHA
        ENDDO
        THPEQ=UTHRS
        WRITE(6,661) N(2,(JMAX+1)/2),KP,UTHRS
      ENDIF
 661  FORMAT(' H+ and He+ reduced; equatorial [H+]=',F10.2
     >   ,1X,' for large Kp=',F4.1,' at UT=',F8.2)
      RETURN
      END 
C:::::::::::::::::::::::::::::: SAVDEN ::::::::::::::::::::::::::::::::
C..... This routine saves all the values to be used after non-convergence in
C..... Odd N solution
      SUBROUTINE SAVDEN(IFLAG)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      COMMON/ODNSAV/ NSAVE(4,401),TISAV(3,401),XISAV(9,401),XVSAV(9,401)
     >  ,USAV(2,401),UESAV(2,401),STRSAV(12,401),RSAV(401)
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)

      DO J=1,JMAX
        DO I=1,4
           NSAVE(I,J)=N(I,J)
        ENDDO
        DO I=1,3
          TISAV(I,J)=TI(I,J)
        ENDDO
        DO I=1,2
           UESAV(I,J)=UE(I,J)
           USAV(I,J)=U(I,J)
        ENDDO
        DO I=1,12
           STRSAV(I,J)=STR(I,J)
        ENDDO
        RSAV(J)=RCON(J)
        XISAV(1,J)=XIONN(1,J)
        XVSAV(1,J)=XIONV(1,J)
      ENDDO
      RETURN
      END
C:::::::::::::::::::::::::::::: SAVDEN ::::::::::::::::::::::::::::::::
C..... This routine recovers all the values to be used after non-convergence in
C..... Odd N solution
      SUBROUTINE RECDEN(IFLAG)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      COMMON/ODNSAV/ NSAVE(4,401),TISAV(3,401),XISAV(9,401),XVSAV(9,401)
     >  ,USAV(2,401),UESAV(2,401),STRSAV(12,401),RSAV(401)
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)

      IFLAG=-9       !... flag indicates retrieving stored values
      DO J=1,JMAX
        DO I=1,4
           N(I,J)=NSAVE(I,J)
        ENDDO
        DO I=1,3
          TI(I,J)=TISAV(I,J)
        ENDDO
        DO I=1,2
           UE(I,J)=UESAV(I,J)
           U(I,J)=USAV(I,J)
        ENDDO
        DO I=1,12
           STR(I,J)=STRSAV(I,J)
        ENDDO
        RCON(J)=RSAV(J)
        XIONN(1,J)=XISAV(1,J)
        XIONV(1,J)=XVSAV(1,J)
      ENDDO
      RETURN
      END
C::::::::::::::::::::::::: EDEN170 ::::::::::::::::::::
      SUBROUTINE EDEN170(JMAX,     !.. # points on grid
     >                   Z,        !.. Altitude array
     >                   N,        !.. Ion density array
     >                   J170,     !.. index of 170 km alt in Z
     >                   N170KMN,  !.. [e] at 170 km in north
     >                   N170KMS)  !.. [e] at 170 km in south
      IMPLICIT NONE
      INTEGER J,JMAX,J170,J170C
      DOUBLE PRECISION Z(401),N(4,401),N170KMN,N170KMS
      J170=2 
      DO J=2,JMAX/2
        IF(170.LT.(Z(J+1)+Z(J))/2.AND.170.GE.(Z(J)+Z(J-1))/2) J170=J
      ENDDO
      N170KMN=N(1,J170)+N(2,J170)+N(3,J170)+N(4,J170)
      J170C=JMAX-J170+1
      N170KMS=N(1,J170C)+N(2,J170C)+N(3,J170C)+N(4,J170C)
      RETURN
      END
C::::::::::::::::::::::::: EDEN170 ::::::::::::::::::::
C.. This routine calculates the sum of the squares of the differences between the
C.. measured and modeled hmF2
      SUBROUTINE CAL_STATS(SEC,DT,HMF2N,DHMF2N,HMF2S,DHMF2S)
      IMPLICIT NONE
      INTEGER JTI
      REAL SEC
      DOUBLE PRECISION HMF2N,DHMF2N,HMF2S,DHMF2S,SUMSQN,SUMSQS,NUMSQ,DT,
     >  SUMN,SUMS
      DATA JTI/0/

      IF(JTI.GT.-1.1111) RETURN  !...

      JTI=JTI+1
      !.. zero out every 24 hours to see problem data
      IF(JTI.EQ.1.OR.AMOD(SEC,86400.0).LT.DT+1) THEN
       IF(NUMSQ.GT.0) WRITE(6,441) 
     >   NINT(NUMSQ),DSQRT(SUMSQN/NUMSQ),DSQRT(SUMSQS/NUMSQ),
     >   SUMN/NUMSQ,SUMS/NUMSQ
 441    FORMAT(' #=',I7,' SIGMAN=',1PE10.2,' SIGMAS=',1PE10.2,
     >   ' SUMN=',1PE10.2,' SUMS=',1PE10.2)
        SUMSQN=0.0
        SUMSQS=0.0
        NUMSQ=0.0
        SUMN=0.0
        SUMS=0.0
      ENDIF

      NUMSQ=NUMSQ+1
      SUMSQN=SUMSQN+(DHMF2N-HMF2N)**2
      SUMSQS=SUMSQS+(DHMF2S-HMF2S)**2
      SUMN=SUMN+(DHMF2N-HMF2N)
      SUMS=SUMS+(DHMF2S-HMF2S)

      RETURN
      END
C:::::::::::::::::::::::::: ADJDTMAX :::::::::::::::::::::
C..... This subroutine is used with ExB drift to shorten time step near the end
C..... of the run so that the run ends as close as possible to specified time
C..... Written by P. Richards December 2008
      SUBROUTINE ADJDTMAX(ITAU,     !.. Total run time
     >                      DT,     !.. Current time step
     >                   TSTOP,     !.. MLength of run
     >                   TWIDT,     !.. Twilight time step
     >                   DTMAX)     !.. Maximum time step
      IMPLICIT NONE
      INTEGER ITAU
      DOUBLE PRECISION DT,TSTOP,TWIDT,DTMAX
      !.. Don't allow max time step to be too big for ExB
      IF(DTMAX.GT.3000) THEN
         WRITE(6,*) ' ExB drift: Max time step reset to 3000 secs'
         DTMAX=3000
      ENDIF
      !.. Shorten time step near end of the run to end closer to end L
      IF((REAL(ITAU)+DT)/3600.0.GT.TSTOP-0.5) THEN
        IF(DTMAX.GT.300) DTMAX=300
        IF((REAL(ITAU)+DT)/3600.0.GT.TSTOP-0.1.AND.DTMAX.GT.60) DTMAX=60
        TWIDT=DTMAX   !.. make sure the results are stored in PRINT file
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C...... This subroutine calculates the limiting H+ flux in the Northern Hemisphere
C...... using the method described by Richards and Torr JGR 1985 page 5261
C...... Updated July 2017
      SUBROUTINE FLXLIMIT(JR,  !.. Index of limiting flux starting altitude. See below
     >                    JQ,  !.. Index of equatorial grid point
     >                    JK,  !.. Index of NmF2 grid point
     >                     Z,  !.. Altitude grid (km)
     >                     N,  !.. Ion densitie array O+, H+, tot+, He+ (cm-3)
     >                    TI,  !.. Ion and electron temperatures (K)
     >                    BM,  !.. Magnetic field strength
     >                    GR,  !.. Gravity array (cm-2s-1)
     >                    ON,  !.. Atomic oxygen density array (cm-3)
     >                    HN,  !.. Atomic hydrogen density array (cm-3)
     >                    TN,  !.. Neutral temperature array (K)
     >                    SL,  !.. Distance along field line (cm)
     >                    JZ,  !.. OUT: Index of ZTRAN
     >                 FYLIM,  !.. OUT: Limiting flux (cm-2 s-1)
     >                HPLOSS,  !.. OUT: H+ total loss above escape altitude (cm-2)
     >                    HS,  !.. OUT: O+ scale height (cm)
     >                 ZTRAN)  !.. OUT: Altitude of transition to escape flux (km)
      INTEGER JR,JQ,JK
      DOUBLE PRECISION RTS(199),ON(401),HN(401),TN(401),BM(401),SL(401)
     > ,Z(401),N(4,401),TI(3,401),GR(2,401),H1,H2,H0,ACON,ZTRAN,HS,
     > FYLIM,HPLOSS
     
      !.. JR is the grid point corresponding to the reference altitude (Zr) 250 to
      !.. 300km below the expected transition altitude. FLIP uses
      !.. Zr = 200.0+200.0*(F107+F107A)/2.0/70.

      !.. Determine the H+ flux by simple method
      H1=-51.6*(TI(2,JR)+TI(3,JR))/GR(2,JR)  !.. O+ scale height
      H2=-51.6*TN(JR)/GR(2,JR)               !.. O scale height
      H0=1.0E5*H1*H2/(H1-H2)
      
      !.. H+ + O and O+ + H reaction rates are from Fox and Sung JGR 2001, page 21,305
      CALL RATS(JR,TI(3,JR),TI(2,JR),TN(JR),RTS)  !.. get reaction rates

      !.. Calculate the transition altitude
      ACON=9E7*TI(2,JR)**2.5/RTS(1)/H0**2
      ZTRAN=Z(JR)-DLOG(ACON/N(1,JR)/ON(JR))*H1*H2/(H1+H2)

      !.. Find the index of the transition altitude index JZ
      JZ=JK   !.. starting at hmF2 altitude and ascending
      DO J=1,JQ
        IF(J.GT.JK.AND.Z(J).LE.ZTRAN) JZ=J
      ENDDO

C       !.. Alternative for finding transition altitude. Use the peak H+ density.
C       It is lower than ZTRAN and overestimates the limiting flux by about 25%
C      JZ=JQ   !.. Begin at equator and go down the field line towards hmF2
C      DO J=JQ-1,JK,-1
C        !.. Check for peak density altitude
C        IF(N(2,J).GT.N(2,J-1).AND.N(2,J).GT.N(2,J+1))
C     >    JZ=J
C      ENDDO

      IF(ABS(Z(JZ)-Z(JR)).GT.500.0) JZ=JR  !.. error capture
      
      HS=-51.6E5*(TI(2,JZ)+TI(3,JZ))/GR(2,JZ)      !.. O+ scale height
      !.. Limiting flux normalized to the first grid point cross section
      !.. For comparison, the FLIP flux is also normalized to the same 
      !.. altitude for output
      FYLIM=RTS(2)*HN(JZ)*N(1,JZ)*HS*BM(1)/BM(JZ)

      !.. Integrate the H+ total loss above ZTRAN. Should be small
      HPLOSS=0.0D0
      DO J=JZ,JQ
        HPLOSS=HPLOSS+RTS(1)*N(2,J)*ON(J)*(SL(J+1)-SL(J))
      ENDDO
      RETURN
      END
C:::::::::::::::::::::::::::::: FLIP_MODS ::::::::::::::::::::::::::::::::::::::::
C......This routine records changes made since May 2017
      SUBROUTINE FLIP_MODS
      IMPLICIT NONE
      WRITE(6,'(//A)') ' ====================================='
      WRITE(6,'(/A,//3X,A/)') ' ..... FLIP model Modifications .....',
     >  'Search for the date in the FLIP model to find the modification'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2017-05-20: Modified routine ALTWIN to reduce wind speeds'
     > ,' when the input hmF2 is below normal F2 altitudes (<~215km).'
     > ,' This is mainly a problem when there is no F2 peak, which'
     > ,' is known as the G condition on ionograms. It is'
     > ,' fairly common around midday in summer at solar minimum,'
     > ,' when O+ loss rates are high.'
     > ,' FLIP uses a running mean for missing hmF2, but this may not'
     > ,' represent the actual situation. As a result, the hmF2 winds'
     > ,' can become excessive and erratic. The modification reduces'
     > ,' the number and size of excessive winds.'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2017-07-08: Redid the Richards and Torr JGR 1985 limiting'
     >  ,' flux.  Results are in the end columns of the *Heat*.txt file'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2017-09-08: Redid solar eclipse routine (ECFUN) to split the'
     >  ,' irradiance into chromosphere responsible for photoinization'
     >  ,' and corona responsible for photoelectrons. Subsequently'
     >  ,' modified to vary 30.4 nm as chromosphere emission'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2018-02-27: Odd N and N2(v) solns switched off for L < 1.07.'
     >  ,' A warning is provided for ExB drift for L < 1.15'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2018-03-16: RATCHK routine modified to test reaction rates'
     >  ,' at 3 different temperatures.'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2018-10-24: Fixed a bug when using a file for EUV.'
     >  ,' Photoionization was not calculated in print phase.'     
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2018-10-30: Updated IRI95 routine for when a requested year'
     >  ,' is greater than the end year in the IG-RZ index file.'    
     >  ,' Now, a message explains how to update the file.'    
      WRITE(6,'(/A,A)') 
     >  ' 2018-11-19: Further updated limiting flux calculation and'
     >  ,' its output to the *HEAT*.txt file.'    
      WRITE(6,'(/A,A)') 
     >  ' 2018-12-03: Made the Kp output real in the *LOG* file.'
     >  ,' Also made it more accurate by using a lookup table'    
      WRITE(6,'(/A,A/2X,A)') 
     >  ' 2020-07-23: Fixed a problem in CALVN2 with N2v factor just'
     >  ,' below the field line grid apex point at low latitudes.'
     >  ,' Also modified vertical grid spacing in VERGRD.'
      WRITE(6,'(/A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2020-07-23: Created routine TEST_WIND_FILES to test wind/hmF2'
     >  ,' file for winds or hmF2. It is used at low latitudes where'
     >  ,' it is not possible to get winds from hmF2. It allows'
     >  ,' winds to be input from a file'
     >  ,'2020-07-29: North winds are positive in both files now.'
      WRITE(6,'(/A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2020-07-29: Modified routine XION in MINORA.FOR to'
     >  ,' chemically calculate He+ when IHEP=0, and N+ when' 
     >  ,' INPLUS=0 in FLIPRUN.DDD'
      WRITE(6,'(/A,A,A)') 
     >  ' 2020-08-02: Fixed a problem in calculating the index of hmF2'
     >  ,' when there is no peak above 180 km in routine DIAGPR.'
      WRITE(6,'(/A,A,/2X,A)') 
     >  ' 2020-08-02: Modified routine CEMOX to output 6300A and 7755A'
     >  ,' emissions for Bob Kerr. TWI-EMIS was updated to accomodate'
     >  ,' the new FLIP-Print output'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-08-30: Changes for getting output above a near equatorial'
     >  ,'station. Multiple runs are made with different L-values'
     >  ,'to produce an electron density profile. Output is on unit 43.'
     >  ,'It is turned on by specifying magnetic latitude with'
     >  ,'the 8th parameter in FLIPRUN.DDD, which was always zero.'
     >  ,'For backward compatibility, if it is < 0.001, no output.'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-09-10: Improved error messages for inputting winds' 
     >  ,'at high and equatorial latitudes'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-09-27: Added HEATFS and TOPTENSV CAL_EHEAT to smooth out'
     >  ,' big jumps in Te data. Also since the Te algorithm needs a Te'
     >  ,' difference to work, a small percentage is added to TOPTEN'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-10-28: Fixed the Ti/Tn bug in Rates.for'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-11-06: Took care of TIMEOF < TIMEON in FLIPRUN.DDD. By'
     >  ,'FLIP adding 24 to TIMEOF for the run.'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2021-02-02: The reduction in plasmaspheric H+ and He+ after'
     >  ,'high Kp was changed to Chappell Rev. Geophys 1972 figure 4.'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)')
     >  '2021-07-15: Implemented Apex coordinates. FLIP may need'
     >  ,'updating after 2025 but only if north and south geographic'
     >  ,'latitude and longitude are too different from the IGRF'
     >  ,'values. To update, replace the Apex.f90 file in the folder'
     >  ,'FLIP with the latest version from the Web and recompile.'
     >  ,'https://github.com/NCAR/apex_fortran/blob/master/apex.f90'
     >  ,'https://omniweb.gsfc.nasa.gov/vitmo/cgm.html'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-05-16: Modified the hmF2 proxy wind algorithm to allow'
     >  ,'hmF2 to go down to 195 km based on Kharkiv December 2019'
     >  ,'radar data showing hmF2 at 200 km'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-05-29: FLIPRUN.DDD modified to allow different H density'
     >  ,'multipliers for the Northern and Southern hemispheres.'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-08-06: Fixed a problem with the field line grid for'
     >  ,'low L values when the field line does not reach 250 km.'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-08-15: Modification to allow Hydrogen densities to be'
     >  ,'read in from the hmF2/wind data files'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-10-15: RSLPSD and RSTEMC_ExB modified to output Te heat'
     >  ,'equation details for Bill Peterson'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-11-01: Modification to print file headers in most files'
     >  ,'every 24 hours at midnight LT'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2023-01-09: FLIP run & print files modified to write altitude'
     >  ,'profiles of HWM14 (zon. & mer.) winds & FLIP magnetic'
     >  ,' meridional winds. Also added winds at a specified altitude'
     >  ,'to the UT LT SZA ... output line at the FLIP-print stage'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2023-02-06: Modified run & print files to process gravity'
     >  ,'wave data from a file'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2023-02-16: Changed UT output from F8.2 to F8.3'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2023-08-19: Modified to add Gravity Waves to winds input'
     >  ,'from a file'
      RETURN
      END